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We use large scale quantum Monte Carlo simulations to study an extended version of the canonical 
Shastry-Sutherland model - including additional interaction and exchange anisotropy - over a wide 
range of Hamiltonian parameters and applied magnetic field. The model is appropriate for describing 
the low energy properties of some members of the rare earth tetraborides. Working in the limit of 
large Ising-like exchange anisotropy, we demonstrate the stabilization of columnar antiferromagnetic 
(CAFM) order in the ground state at zero field and an extended magnetization plateau at m/nis — 
1/2 in the presence of an applied longitudinal magnetic field - qualitatively similar to experimentally 
observed low-temperature phases in ErB4. Our results show that for an optimal range of exchange 
parameters, a spin supersolid ground state is realized over a finite range of applied field between the 
CAFM phase and the magnetization plateau. The full momentum dependence of the longitudinal 
and transverse components of the static structure factor is calculated in the spin supersolid phase to 
demonstrate the simultaneous existence of diagonal and off-diagonal long-range order. Our results 
will provide crucial guidance in designing further experiments to search for the elusive spin supersolid 
phase in ErB4. 
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Ever since Penrose and Onsager [T] speculated on the 
possible coexistence of diagonal and ofF-diagonal long 
range order, supersolid order has been of broad interest 
within the physics community. Although the initial in- 
spiration for this exotic type of ordering originated from 
the consideration of the remarkable properties of solid 
helium, the realization of a supersolid phase of helium 
remains controversial. [3S] On the other hand, theoretical 
studies of several models of lattice bosons with competing 
interactions have conclusively established the presence of 
supersolid phases over extended parameter regimes. 
In this case, the discreteness of the lattice simplifies the 
process of forming a supersolid, and the bosonic mod- 
els can potentially be experimentally realized with cold 
atoms in optical lattices. Concurrently, the realization 
of BEC of magnons and other novel bosonic phases in 
quantum magnets provides an alternative route to this 
elusive state of matter, and the spin analog of the super- 
solid phase was introduced. 0[S]. Yet to date, no experi- 
mental system has been found that displays unambiguous 
signs of supersolid ordering. Hence, the identification of 
realistic systems with supersolid order is a valuable step 
towards the ultimate realization of this elusive phase. 

A key element in the stabilization of a supersolid 
ground state in many models is geometric frustration, 
which makes frustrated quantum magnets the natu- 
ral place to look for the spin supersolid phase. The 
rare-earth tetraborides (RB4) are a promising class of 
quantum magnets for studying the effects of geometric 
frustration in interacting spin systems. Consisting of 
weakly coupled layers of magnetic moment carrying R'^^" 
ions arranged on the Shastry-Sutherland lattice (SSL), 
their magnetic properties are well described by appro- 
priate generalizations of the Shastry Sutherland model 



(SSM). [5] Several members of this family have been ob- 
served to exhibit magnetization plateaus at low temper- 
atures. For example, both ErB4 and TmB4 exhibit an 
extended magnetization plateau at 1/2 the saturation 
magnetization. [lOl [11] Despite this similarity, however, 
additional fractional magnetization plateaus have been 
observed in TmB4 that are not seen in ErB4. Further, 
the zero-field ground state magnetic structure of these 
two RB4 differ significantly. Recent analysis using neu- 
tron scattering experiments has determined the low tem- 
perature magnetic structure of TmB4 to be collinear and 
antiferromagnetic (AFM) with a Q = (tt, tt) ordering of 
the local moments, [HI [T3] i.e. the Ising analog of the 
Neel state. In contrast, it is well-documented that the 
magnetic order of the zero-field ground state in ErB4, 
while collinear, is antiferromagnetic only in one direc- 
tion, i.e. a columnar ordering with either Q = (7r,0) or 
Q = (0,7r). [inilll] 

Previous studies have shown that the low tempera- 
ture magnetic properties of TmB4 can be explained by 
an effective low energy model that is obtained by ex- 
tending the canonical SSM to incorporate ferromagnetic 
(FM) exchange with Ising-like anisotropy and additional 
long range interaction. [THl US] The additional interac- 
tion is instrumental in stabilizing an extended magneti- 
zation plateau at m/rris = 1/2, while at the same time 
eliminating the magnetization plateau at m/rris = 1/3 
that is ubiquitous to the canonical SSM in the Ising 
limit. [T71[TB| Similar generalizations are expected to cap- 
ture the low-energy properties of other members of the 
RB4 family. Thus, it is useful to extend the results of 
this effective model to include magnetization processes 
that begin from zero-field ground states other than the 
Neel-like antiferromagnetic state considered so far. 



In this Letter, we determine the parameter regime of 
the effective model for RB4 for which a columnar an- 
tiferromagnetic ordering is stabilized in the absence of 
an external field. We present large scale quantum Monte 
Carlo (QMC) simulations of the magnetization process in 
this regime and demonstrate the appearance of a field- 
induced spin supersolid phase at magnetizations below 
half saturation. Arguing that similar factors may drive 
the stabilization of columnar order in ErB4, we identify 
ErB4 as a promising candidate for the observation of su- 
persolid order. 

Many RB4 are well-suited to numeric investigation of 
their low-temperature magnetic properties. The high to- 
tal angular momentum J of the ions, combined with 
a large easy-axis single-ion anisotropy, allows us to con- 
struct a simple effective model of the low energy states. 
Briefly, the large easy-axis single- ion anisotropy, Z?, splits 
the local Hilbert space at each site into degenerate dou- 
blets - the successive doublets are separated by energy 
gaps ~ Z?. To a first approximation, then, the low-energy 
magnetic interactions are captured by an Ising model 
comprised of the maximal- J doublet. Higher order pro- 
cesses lead to weak FM (AFM) transverse exchange in- 
teractions between neighboring sites for integer (half-odd 
integer) values of J, and thus the XXZ model with Ising- 
like exchange anisotropy becomes the preferred model to 
study the low-energy magnetic properties of ErB4 and 
TmB4, both of which possess a large easy- axis single- ion 
anisotropy. |29j In this work, we restrict ourselves to the 
case of FM transverse exchange interactions, where the 
model is free of the negative sign problem and thereby 
directly amenable to QMC simulation. 

The above model can be described by the Hamiltonian 
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where the summation is over all bond types a of strength 
Jq, and their associated bonds {ij)a acting on sites i and 
j of the extended SSL depicted in Fig. [TJa). Here, in 
addition to the canonical interactions Ji and J2 of the 
SSM, we have included an additional interaction J3 along 
the diagonals of the "empty" plaquettes (i.e. those miss- 
ing the J2 interaction). We work in the experimentally 
relevant limit of strong Ising-like exchange anisotropy, 
A <C 1, and consider the effect of an external magnetic 
field hz- Since the crystal structure of these materials is 
such that the bonds of the canonical SSL are expected 
to be of roughly equal strength, in what follows we set 
Ji = J2 = I7 and all parameters are given in these units. 

To characterize the different phases of H, longitudinal 
and transverse components of the static structure factor 




FIG. 1: (Color online) Ising ground state structure of the half 
plateau state for antiferromagnetic J3. Black and white circles 
represent down and up spins, respectively, while blue outlines 
indicate spins that can be flipped with minimal energy cost. 
The bonds Ji, J2, and J3 (see text) are illustrated in (a) as 
solid blue, thick green, and dotted red lines, respectively. In 
(b) the hopping processes in the spin supersolid phase are 
illustrated (J3 bonds are omitted for clarity). 



are defined as 



S'^^(fc) = ^^(S'fS'De*^''--^^) 



(2) 



In the basis, these serve as measures of off-diagonal 
and diagonal order, respectively. As the uniform magne- 
tization per site is simply ruz = '^j{Sj)/N, here we note 
that 5*^^(0) = N{m1). Further, we define staggered and 
columnar magnetizations as {m^) = S^^{t:,'k) /N and 
{ml) = [S'^^(0,7r) + S'^^(7r,0)] /TV, respectively. 

QMC simulations are performed using the stochas- 
tic series expansion algorithm with directed loop up- 
dates. The correlations {S^ SJ + S~ Sp and 
{S^S^) can be straightforwardly computed within this 
method [211 El]. A useful observable in characteriz- 
ing the ground state phases in numerical simulations 
is the spin stiffness, ps, defined as the response to a 
twist in the boundary conditions. [53] In simulations 
that sample multiple winding number sectors, the eval- 
uation of the stiffness simplifies to calculating the wind- 
ing number of the world lines, leading to the expression 
Ps = {wl + w. 



/4/3 in two dimensions, [24j where 



and Wy are the total winding in the x and y directions. 

Let us begin by considering the ground-state phase di- 
agram of % in the Ising limit, A = 0. In this limit, we can 
construct the various possible ground states by hand. By 
comparing their respective energies, the phase diagram 
in the hz — J3 plane is determined (see solid black lines 
and labeled phases in Fig. [2]) . This will provide a back- 
ground of classically ordered magnetic states upon which 
we shall study the effect of quantum fluctuations. At zero 
field the ground state is an Ising antiferromagnet with ei- 
ther staggered (Q = (7r,7r)) or columnar (Q = (7r,0) or 
Q — (0, 7r)) order. We will refer to these states as SAFM 
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and CAFM, respectively. For small J3 there is a field- 
induced magnetic phase transition to a 1/3 plateau (TP) 
with rUz = TOs/3. With increasing field, this state further 
evolves into a 1/2 plateau (HP) with = ms/2, except 
at J3 = 0, where there is a direct transition from the TP 
to the fully polarized (FP) state (m^ = nis)- [21] At large 
absolute values of J3 the magnetization process skips the 
TP, with a direct transition from the zero-field ground 
states to the HP. Notice that there are two separate HP 
states with distinct structures for positive and negative 
J3 values. |26j| For J3 < the HP consists of diagonal 
stripes with weak staggered order and no columnar order, 
while for J3 > the HP is depicted in Fig. [T] and consists 
of alternating FM and AFM stripes. The freedom of the 
AFM stripes to shift phases with respect to one another 
leads to a large degeneracy that eliminates staggered or- 
der, though weak columnar order remains. Thus, the 
two HP states mirror the order of the underlying zero- 
field ground states (SAFM and CAFM, respectively). In 
the remainder of this work, we refer exclusively to the 
HP state with columnar ordering. 

Turning to the case of finite A = 0.10, QMC simula- 
tions are performed to determine the ground state quan- 
tum phase diagram (see Fig. [5]) . The numerically deter- 
mined phase boundaries are marked by data points, while 
three emergent phases with transverse ferromagnetic or- 
der (i.e. superfluid) are delineated by dashed lines. These 
modifications are primarily for antiferromagnetic values 
of J3. For ferromagnetic J3 practically no change is ob- 
served with the addition of small transverse terms. This 
is because the magnetic phase transitions in this region 
are strongly first order, and thus larger quantum fiuc- 
tuations are required to alter the magnetization process. 
The first emergent phase (red squares, Fig.|2| begins from 
weakly ferromagnetic J3 values and extends throughout 
the antiferromagnetic J3 region of the phase diagram, 
where a superfluid (SF) phase opens up between the HP 
and FP phases. A similar phase opens up between the 
CAFM and HP phases for larger values of J3. This sec- 
ond emergent phase (blue diamonds) turns out to be a 
spin supersolid (SS) phase, possessing both superfluid 
and columnar antiferromagnetic order. A third emergent 
phase (green triangles) with superfluid ordering occurs at 
the border of the SAFM, CAFM, and TP phases. 

The phase boundary between the SF and FP states can 
easily be determined by considering the energy required 
to flip a single spin within the FP phase. This energy is 
zero at the boundary, which is given by the relation hs — 
(5/2 -I- J3) (1 -I- A). A similar construction can be used 
to determine the boundary of the SS and HP phases and 
yields he ^2A + (1/2 + J3) (1 + A). We have employed 
these relations in drawing the SF-FP and SS-HP phase 
boundaries shown in Fig. [2j Note, however, that the 
CAFM-SS and HP-SF phase boundaries are drawn as 
fits to the QMC data, and are merely guides to the eye. 

Focusing on the intermediate J3 > regime, in Fig. [3] 
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FIG. 2: (Color online) Phase diagram in the Ising limit (solid 
black lines) along with modifications for the XXZ model 
at Ising-like anisotropy A = 0.10 (colorful dashed lines). 
Ising phase boundaries have been calculated exactly through 
ground state energy evaluation. For the XXZ model, QMC 
boundaries are indicated by data points. 



we show the magnetization process for J3 = 1.25 and 
A ~ 0.10. Here, the zero-field ground state is the CAFM 
phase with near its maximal value of 0.25. As the 
magnetic field is increased, there is a field-induced first- 
order phase transition from the CAFM phase to the SS 
phase, as evidenced by discontinuous jumps in the mag- 
netic observables (see Fig. [3|. In the SS phase, ml is 
reduced, yet remains finite, and becomes non-zero, 
even as a non-zero spin stiffness develops. These are 
signs of supersolid order, which we carefully consider in 
the next paragraph. After a continuous transition to the 
HP state, finite size scaling shows approaching zero, 
despite the disappearance of superfiuid order. This is a 
consequence of the large degeneracy in the HP phase that 
is not present in the SS phase. Next, a first-order phase 
transition (with accompanying jumps in the magnetic ob- 
servables) eliminates the remaining diagonal magnetic or- 
der (to^), as the HP gives way to the SF phase with only 
off-diagonal ordering. Finally, there is a continuous phase 
transition from the SF to the FP state. 

A supersolid is defined by the simultaneous occurrence 
of diagonal and off-diagonal long-range order. To confirm 
that the present spin supersolid phase satisfies this prop- 
erty, we calculate the full static structure factors of the 
transverse and longitudinal degrees of spin freedom (5'' 
and S^^, respectively). As seen in Fig. |4] it is clear that 
both structure factors possess ordering. The transverse 
ordering is primarily at A; = (0, 0) while the longitudi- 
nal ordering is primarily at fc = (0, vr) and k — (7r,0). 
The secondary peak in at k = (7i',7r) is due to the 
off-diagonal order being modulated by the diagonal or- 
der. In contrast, the secondary peak in S^^ is located 
at the origin and simply reflects the fact that the uni- 
form magnetization is non-zero. Finally, a tertiary peak 
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FIG. 3: (Color online) Behavior of the columnar (m^), stag- 
gered (ml), and uniform (mz) magnetizations, along with the 
spin stiffness (ps), for the complete magnetization process at 
A = 0.10, J3 — 1.25, and inverse temperature P = 64. The 
sequence of phases with increasing magnetic field is columnar, 
supersolid, half plateau, superfluid, and fully polarized. 



in S'^^ develops spontaneously at k — (tt.tt) in the SS 
phase through a quantum order-by-disorder process that 
we discuss below. 

The mechanism for superflow in the SS phase can fur- 
ther explain the observed features of the transverse struc- 
ture factor. Since the SS phase occurs at magnetizations 
just below half saturation, the mechanism for supersoUd 
formation can be determined by considering the hopping 
processes available when the HP state is doped with down 
spins. The lowest potential energy doping sites in the HP 
phase for J3 > 1 are depicted by white circles with blue 
outlines in Fig. [ij It is clear to see that down spins at 
these sites can easily delocalize by first order processes 
when the phase of AFM stripes is conserved. In con- 
trast, along the domain wall type structure of out-of- 
phase stripes, only second order hopping processes are 
allowed. These two hopping processes are depicted in 
Fig. Qb). Interestingly, since the domain walls impart 
no potential energy gain, the total energy will be lower 
with the first order hopping processes of phase-matched 
stripes. Thus, the degeneracy of the HP state is partially 
lifted by quantum fluctuations in the SS phase, leading to 
a higher degree of diagonal order. This effect is expected 
to grow stronger with increasing J3 values, and is an ex- 
ample of a quantum "order by disorder" process. pTH^U] 

The model Hamiltonian H can be mapped onto a sys- 
tem of hardcore bosons with nearest-neighbor hopping 
i = — I'/aAj /2 and repulsion V = Ja- For such a system 
on the triangular lattice it has been shown [^-|^ that a 
supersolid phase is stabilized for weak hopping t <^ V 
or A ^ 1. Similarly, a square lattice model of hardcore 
bosons with next-nearest-neighbor repulsion (similar to 
our J2 and J3 terms) is known to possess both colum- 





FIG. 4: (Color online) Transverse and longitudinal structure 
factors in the spin supersolid phase, with S^^{k) and S^^{k) 
as defined in Eq. [2] Data are taken from a 24x24 lattice at 
inverse temperature /3 — 100. The free parameters are set to 
A — 0.10, hz = 1.75 and J3 — 1.25 to assure that we are well 
within the spin supersolid phase as outlined in Fig. [2] 



nar order and supersolid phases (see |31j and references 
therein). While no such supersolid state has been found 
for hardcore bosons on the canonical SSL, [35] we have 
shown that the addition of antiferromagnetic J3 bonds 
allows for the formation of a supersolid phase on the ex- 
tended SSL. 

Let us now turn to the experimentally determined 
magnetization process of the rare-earth tetraboride ErB4. 
Our calculations have shown that it is possible to obtain 
the observed magnetization process starting from a zero- 
field columnar antiferromagnetic state. The theoretically 
calculated magnetization process compares well not only 
qualitatively, but also quantitatively, as reflected by the 
extent of the HP state {Ahi/2/h, - (3.5-2)/4 or 38% in 
our model vs. A/ii/a/Zi^ ~ (4-2)/5 or 40% in Ref. [lO]). 
Additionally, this magnetization process may contain a 
spin supersolid phase for magnetizations below half sat- 
uration, as determined by our simulations of an effective 
model with FM transverse exchange interactions. While 
the total angular momentum J — 15/2 for ErB4 actu- 
ally leads to AFM transverse exchange interactions, we 
believe that a similar mechanism could lead to a spin su- 
persolid phase with AFM transverse order. Interestingly, 
such an ordering in the effective transverse spin degrees 
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of freedom would correspond to antiferronematic order 
of the total angular momentum. [33' 
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